source("include.R")
❤ Loaded agrimonia dataset. Available as df_agri.
❤ Converted df_agri$Time to Date variable type, year-month-day.
❤ Created DATE_FORMAT variable for date comparisons.
❤ Loaded stations split dataset. Available as df_stat.
Maps plots
library(maps)
Avvertimento: il pacchetto ‘maps’ è stato creato con R versione 4.2.3
library(ggplot2)
Caricamento pacchetto: ‘ggplot2’
Il seguente oggetto è mascherato da ‘package:crayon’:
%+%
library(sf)
Avvertimento: il pacchetto ‘sf’ è stato creato con R versione 4.2.3Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
Idea: plot the values for pm10 each day, like with “bubbles” of
ggplot, so to get 2192 plots. And then combine all them in a single gif,
which shows the evolution over time of those pm10 values.
https://askubuntu.com/questions/648244/how-do-i-create-an-animated-gif-from-still-images-preferably-with-the-command-l
So idea: plot bubbles with size according to the pm10 value. And
maybe if we want use also a different color (like to separate
high/medium/low values).
regioni_italiane <- st_read("italia/gadm40_ITA_1.shp")
Reading layer `gadm40_ITA_1' from data source
`C:\Users\feder\Desktop\Uni magistrale\Bayesian statistics\progetto-bayesian\src\italia\gadm40_ITA_1.shp'
using driver `ESRI Shapefile'
Simple feature collection with 20 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 6.630879 ymin: 35.49292 xmax: 18.52069 ymax: 47.09096
Geodetic CRS: WGS 84
sites <- data.frame(
longitude = unique(df_agri$Longitude),
latitude = unique(df_agri$Latitude))
pad = 0.2*c(-1,1)
mappa_completa <- ggplot() +
geom_sf(data = regioni_italiane, fill = "lightblue", color = "black", size = 0.5)+
coord_sf(xlim = range(sites$longitude)+pad,
ylim = range(sites$latitude)+pad, expand = FALSE)+
# geom_sf(data = world) +
geom_point(data = sites, aes(x = longitude, y = latitude), size = 1, shape = 23, fill = "darkred")
print(mappa_completa)

sites <- data.frame(
longitude = unique(df_agri$Longitude),
latitude = unique(df_agri$Latitude))
for (data in seq.Date(from=as.Date("2016-01-01",DATE_FORMAT),
to=as.Date("2016-01-11",DATE_FORMAT),by="days")){
dt = as.Date(data,origin="1970-01-01")
df_date = df_agri[which(df_agri$Time == dt),]
print(dt)
pad = 0.2*c(-1,1)
mappa_completa <- ggplot() +
geom_sf(data = regioni_italiane, fill = "lightblue", color = "black", size = 0.5)+
coord_sf(xlim = range(sites$longitude)+pad,
ylim = range(sites$latitude)+pad, expand = FALSE)+
# geom_sf(data = world) +
geom_point(data = sites, aes(x = longitude, y = latitude), size = 1, shape = 23, fill = "darkred")+
geom_point(data=df_date,aes(x=Longitude,y=Latitude),size=df_date$AQ_pm10/10,alpha=0.5)+
ggtitle(dt)
print(mappa_completa)
ggsave(paste0("./maps_images/",dt,".jpeg"),plot=mappa_completa)
}
[1] "2016-01-01"
[1] "2016-01-02"
[1] "2016-01-03"
[1] "2016-01-04"
[1] "2016-01-05"
[1] "2016-01-06"
[1] "2016-01-07"
[1] "2016-01-08"
[1] "2016-01-09"
[1] "2016-01-10"
[1] "2016-01-11"











Trend plots
stations = unique(df_agri$IDStations)
cols = colora(6,97,show=1)

for (st in stations[1:4]){
df_st = df_stat[[st]]
data2016 = subset(df_st,Time>="2016-01-01" & Time<="2016-12-31")
data2017 = subset(df_st,Time>="2017-01-01" & Time<="2017-12-31")
data2018 = subset(df_st,Time>="2018-01-01" & Time<="2018-12-31")
data2019 = subset(df_st,Time>="2019-01-01" & Time<="2019-12-31")
data2020 = subset(df_st,Time>="2020-01-01" & Time<="2020-12-31")
data2021 = subset(df_st,Time>="2021-01-01" & Time<="2021-12-31")
plot(as.Date(format(data2016$Time,"%d %m"),"%d %m"),data2016$AQ_pm10 ,type="l",col=cols[1],
xlab="time",ylab="PM_10 values",main=paste0("Station ",st,", all years"),
ylim=range(na.omit(df_st$AQ_pm10)))
lines(as.Date(format(data2017$Time,"%d %m"),"%d %m"),data2017$AQ_pm10,type="l",col=cols[2])
lines(as.Date(format(data2018$Time,"%d %m"),"%d %m"),data2018$AQ_pm10,type="l",col=cols[3])
lines(as.Date(format(data2019$Time,"%d %m"),"%d %m"),data2019$AQ_pm10,type="l",col=cols[4])
lines(as.Date(format(data2020$Time,"%d %m"),"%d %m"),data2020$AQ_pm10,type="l",col=cols[5])
lines(as.Date(format(data2021$Time,"%d %m"),"%d %m"),data2021$AQ_pm10,type="l",col=cols[6])
legend('topright', levels(as.factor(2016:2021)), fill=cols, bty='n', cex = 0.6)
}




cols = colora(10,show=0)
for (year in 2016:2021){
for (i in 1:length(stations)){
st = stations[i]
df_st = df_stat[[st]]
df_st = df_st[df_st$Time>=as.Date(paste0(year,"-01-01"),"%Y-%m-%d") &
df_st$Time<=as.Date(paste0(year,"-12-31"),"%Y-%m-%d"),]
if(st=="1264"){
plot(df_st$Time,df_st$AQ_pm10,type="l",col=cols[i%%10],
ylim=range(na.omit(df_agri$AQ_pm10)),main=paste0(year," all the stations"),
xlab="months",ylab="PM_10 values")
} else{
lines(df_st$Time,df_st$AQ_pm10,type="l",col=cols[i%%10])
}
}
}






LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQpgYGB7cn0NCnNvdXJjZSgiaW5jbHVkZS5SIikNCmBgYA0KDQojIE1hcHMgcGxvdHMNCmBgYHtyfQ0KbGlicmFyeShtYXBzKQ0KbGlicmFyeShnZ3Bsb3QyKQ0KbGlicmFyeShzZikNCmBgYA0KDQpJZGVhOiBwbG90IHRoZSB2YWx1ZXMgZm9yIHBtMTAgZWFjaCBkYXksIGxpa2Ugd2l0aCAiYnViYmxlcyIgb2YgZ2dwbG90LCBzbyB0byBnZXQgMjE5MiBwbG90cy4gQW5kIHRoZW4gY29tYmluZSBhbGwgdGhlbSBpbiBhIHNpbmdsZSBnaWYsIHdoaWNoIHNob3dzIHRoZSBldm9sdXRpb24gb3ZlciB0aW1lIG9mIHRob3NlIHBtMTAgdmFsdWVzLg0KDQpodHRwczovL2Fza3VidW50dS5jb20vcXVlc3Rpb25zLzY0ODI0NC9ob3ctZG8taS1jcmVhdGUtYW4tYW5pbWF0ZWQtZ2lmLWZyb20tc3RpbGwtaW1hZ2VzLXByZWZlcmFibHktd2l0aC10aGUtY29tbWFuZC1sDQoNClNvIGlkZWE6IHBsb3QgYnViYmxlcyB3aXRoIHNpemUgYWNjb3JkaW5nIHRvIHRoZSBwbTEwIHZhbHVlLg0KQW5kIG1heWJlIGlmIHdlIHdhbnQgdXNlIGFsc28gYSBkaWZmZXJlbnQgY29sb3IgKGxpa2UgdG8gc2VwYXJhdGUgaGlnaC9tZWRpdW0vbG93IHZhbHVlcykuDQoNCmBgYHtyfQ0KcmVnaW9uaV9pdGFsaWFuZSA8LSBzdF9yZWFkKCJpdGFsaWEvZ2FkbTQwX0lUQV8xLnNocCIpIA0KDQpzaXRlcyA8LSBkYXRhLmZyYW1lKA0KCWxvbmdpdHVkZSA9IHVuaXF1ZShkZl9hZ3JpJExvbmdpdHVkZSksIA0KCWxhdGl0dWRlID0gdW5pcXVlKGRmX2FncmkkTGF0aXR1ZGUpKQ0KDQpwYWQgPSAwLjIqYygtMSwxKQ0KbWFwcGFfY29tcGxldGEgPC0gZ2dwbG90KCkgKw0KICBnZW9tX3NmKGRhdGEgPSByZWdpb25pX2l0YWxpYW5lLCBmaWxsID0gImxpZ2h0Ymx1ZSIsIGNvbG9yID0gImJsYWNrIiwgc2l6ZSA9IDAuNSkrDQoJICBjb29yZF9zZih4bGltID0gcmFuZ2Uoc2l0ZXMkbG9uZ2l0dWRlKStwYWQsIA0KCSAgCQkgeWxpbSA9IHJhbmdlKHNpdGVzJGxhdGl0dWRlKStwYWQsIGV4cGFuZCA9IEZBTFNFKSsNCiAgIyBnZW9tX3NmKGRhdGEgPSB3b3JsZCkgKw0KICBnZW9tX3BvaW50KGRhdGEgPSBzaXRlcywgYWVzKHggPSBsb25naXR1ZGUsIHkgPSBsYXRpdHVkZSksIHNpemUgPSAxLCBzaGFwZSA9IDIzLCBmaWxsID0gImRhcmtyZWQiKQ0KIA0KcHJpbnQobWFwcGFfY29tcGxldGEpDQpgYGANCg0KDQpgYGB7cn0NCnNpdGVzIDwtIGRhdGEuZnJhbWUoDQoJbG9uZ2l0dWRlID0gdW5pcXVlKGRmX2FncmkkTG9uZ2l0dWRlKSwgDQoJbGF0aXR1ZGUgPSB1bmlxdWUoZGZfYWdyaSRMYXRpdHVkZSkpDQoNCmZvciAoZGF0YSBpbiBzZXEuRGF0ZShmcm9tPWFzLkRhdGUoIjIwMTYtMDEtMDEiLERBVEVfRk9STUFUKSwNCgkJCSAgIHRvPWFzLkRhdGUoIjIwMTYtMDEtMTEiLERBVEVfRk9STUFUKSxieT0iZGF5cyIpKXsNCgkJDQoJZHQgPSBhcy5EYXRlKGRhdGEsb3JpZ2luPSIxOTcwLTAxLTAxIikNCglkZl9kYXRlID0gZGZfYWdyaVt3aGljaChkZl9hZ3JpJFRpbWUgPT0gZHQpLF0NCglwcmludChkdCkNCgkNCglwYWQgPSAwLjIqYygtMSwxKQ0KCW1hcHBhX2NvbXBsZXRhIDwtIGdncGxvdCgpICsNCgkgIGdlb21fc2YoZGF0YSA9IHJlZ2lvbmlfaXRhbGlhbmUsIGZpbGwgPSAibGlnaHRibHVlIiwgY29sb3IgPSAiYmxhY2siLCBzaXplID0gMC41KSsNCgkJICBjb29yZF9zZih4bGltID0gcmFuZ2Uoc2l0ZXMkbG9uZ2l0dWRlKStwYWQsIA0KCQkgIAkJIHlsaW0gPSByYW5nZShzaXRlcyRsYXRpdHVkZSkrcGFkLCBleHBhbmQgPSBGQUxTRSkrDQoJICAjIGdlb21fc2YoZGF0YSA9IHdvcmxkKSArDQoJICBnZW9tX3BvaW50KGRhdGEgPSBzaXRlcywgYWVzKHggPSBsb25naXR1ZGUsIHkgPSBsYXRpdHVkZSksIHNpemUgPSAxLCBzaGFwZSA9IDIzLCBmaWxsID0gImRhcmtyZWQiKSsNCgkJZ2VvbV9wb2ludChkYXRhPWRmX2RhdGUsYWVzKHg9TG9uZ2l0dWRlLHk9TGF0aXR1ZGUpLHNpemU9ZGZfZGF0ZSRBUV9wbTEwLzEwLGFscGhhPTAuNSkrDQoJCWdndGl0bGUoZHQpDQoJDQoJcHJpbnQobWFwcGFfY29tcGxldGEpDQoJZ2dzYXZlKHBhc3RlMCgiLi9tYXBzX2ltYWdlcy8iLGR0LCIuanBlZyIpLHBsb3Q9bWFwcGFfY29tcGxldGEpDQp9DQpgYGANCg0KDQoNCiMgVHJlbmQgcGxvdHMNCmBgYHtyfQ0Kc3RhdGlvbnMgPSB1bmlxdWUoZGZfYWdyaSRJRFN0YXRpb25zKQ0KY29scyA9IGNvbG9yYSg2LDk3LHNob3c9MSkNCg0KZm9yIChzdCBpbiBzdGF0aW9uc1sxOjRdKXsNCglkZl9zdCA9IGRmX3N0YXRbW3N0XV0NCglkYXRhMjAxNiA9IHN1YnNldChkZl9zdCxUaW1lPj0iMjAxNi0wMS0wMSIgJiBUaW1lPD0iMjAxNi0xMi0zMSIpDQoJZGF0YTIwMTcgPSBzdWJzZXQoZGZfc3QsVGltZT49IjIwMTctMDEtMDEiICYgVGltZTw9IjIwMTctMTItMzEiKQ0KCWRhdGEyMDE4ID0gc3Vic2V0KGRmX3N0LFRpbWU+PSIyMDE4LTAxLTAxIiAmIFRpbWU8PSIyMDE4LTEyLTMxIikNCglkYXRhMjAxOSA9IHN1YnNldChkZl9zdCxUaW1lPj0iMjAxOS0wMS0wMSIgJiBUaW1lPD0iMjAxOS0xMi0zMSIpDQoJZGF0YTIwMjAgPSBzdWJzZXQoZGZfc3QsVGltZT49IjIwMjAtMDEtMDEiICYgVGltZTw9IjIwMjAtMTItMzEiKQ0KCWRhdGEyMDIxID0gc3Vic2V0KGRmX3N0LFRpbWU+PSIyMDIxLTAxLTAxIiAmIFRpbWU8PSIyMDIxLTEyLTMxIikNCgkNCgkNCglwbG90KGFzLkRhdGUoZm9ybWF0KGRhdGEyMDE2JFRpbWUsIiVkICVtIiksIiVkICVtIiksZGF0YTIwMTYkQVFfcG0xMCAsdHlwZT0ibCIsY29sPWNvbHNbMV0sDQoJCSB4bGFiPSJ0aW1lIix5bGFiPSJQTV8xMCB2YWx1ZXMiLG1haW49cGFzdGUwKCJTdGF0aW9uICIsc3QsIiwgYWxsIHllYXJzIiksDQoJCSB5bGltPXJhbmdlKG5hLm9taXQoZGZfc3QkQVFfcG0xMCkpKQ0KCWxpbmVzKGFzLkRhdGUoZm9ybWF0KGRhdGEyMDE3JFRpbWUsIiVkICVtIiksIiVkICVtIiksZGF0YTIwMTckQVFfcG0xMCx0eXBlPSJsIixjb2w9Y29sc1syXSkNCglsaW5lcyhhcy5EYXRlKGZvcm1hdChkYXRhMjAxOCRUaW1lLCIlZCAlbSIpLCIlZCAlbSIpLGRhdGEyMDE4JEFRX3BtMTAsdHlwZT0ibCIsY29sPWNvbHNbM10pDQoJbGluZXMoYXMuRGF0ZShmb3JtYXQoZGF0YTIwMTkkVGltZSwiJWQgJW0iKSwiJWQgJW0iKSxkYXRhMjAxOSRBUV9wbTEwLHR5cGU9ImwiLGNvbD1jb2xzWzRdKQ0KCWxpbmVzKGFzLkRhdGUoZm9ybWF0KGRhdGEyMDIwJFRpbWUsIiVkICVtIiksIiVkICVtIiksZGF0YTIwMjAkQVFfcG0xMCx0eXBlPSJsIixjb2w9Y29sc1s1XSkNCglsaW5lcyhhcy5EYXRlKGZvcm1hdChkYXRhMjAyMSRUaW1lLCIlZCAlbSIpLCIlZCAlbSIpLGRhdGEyMDIxJEFRX3BtMTAsdHlwZT0ibCIsY29sPWNvbHNbNl0pDQoJbGVnZW5kKCd0b3ByaWdodCcsIGxldmVscyhhcy5mYWN0b3IoMjAxNjoyMDIxKSksIGZpbGw9Y29scywgYnR5PSduJywgY2V4ID0gMC42KQ0KDQp9DQpgYGANCg0KDQpgYGB7cn0NCmNvbHMgPSBjb2xvcmEoMTAsc2hvdz0wKQ0KZm9yICh5ZWFyIGluIDIwMTY6MjAyMSl7DQoJZm9yIChpIGluIDE6bGVuZ3RoKHN0YXRpb25zKSl7DQoJCXN0ID0gc3RhdGlvbnNbaV0NCgkJZGZfc3QgPSBkZl9zdGF0W1tzdF1dDQoJCWRmX3N0ID0gZGZfc3RbZGZfc3QkVGltZT49YXMuRGF0ZShwYXN0ZTAoeWVhciwiLTAxLTAxIiksIiVZLSVtLSVkIikgJg0KCQkJCQkgIGRmX3N0JFRpbWU8PWFzLkRhdGUocGFzdGUwKHllYXIsIi0xMi0zMSIpLCIlWS0lbS0lZCIpLF0NCgkJaWYoc3Q9PSIxMjY0Iil7DQoJCQlwbG90KGRmX3N0JFRpbWUsZGZfc3QkQVFfcG0xMCx0eXBlPSJsIixjb2w9Y29sc1tpJSUxMF0sDQoJCQkJIHlsaW09cmFuZ2UobmEub21pdChkZl9hZ3JpJEFRX3BtMTApKSxtYWluPXBhc3RlMCh5ZWFyLCIgYWxsIHRoZSBzdGF0aW9ucyIpLA0KCQkJCSB4bGFiPSJtb250aHMiLHlsYWI9IlBNXzEwIHZhbHVlcyIpDQoJCX0gZWxzZXsNCgkJCWxpbmVzKGRmX3N0JFRpbWUsZGZfc3QkQVFfcG0xMCx0eXBlPSJsIixjb2w9Y29sc1tpJSUxMF0pDQoJCX0NCgl9DQp9DQpgYGANCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg==